install.packages(c('readxl', 'sf', 'terra', 'leaflet', 'geodata', 'wbstats',
'rnaturalearth', 'osmextract', 'osrm', 'plotly'))GEO UNIV’R Tunisie 2024
LUN3 - Acquisition de données géographiques et visualisations de base
Note
Ce support est inspiré en de nombreux points du manuel Géomatique avec R (Giraud et Pecout 2024). Et on les en remercie chaleureusement 🙏.
Ce document montre les solutions techniques pour importer, explorer, manipuler, visualiser (simplement) et exporter des données géographiques avec R :
- En partant de rien.
- Avec des tableaux de données comprenant une dimension géographique (identifiant géographique).
- Avec des couches d’information géographique (vecteur et raster).
- En utilisant des solutions des packages interfaçant des API (Application Programming Interface).
- En mobilisant des packages interfaçant OpenStreetMap (OSM).
Les exemples proposés sont tantôt adaptés aux contexte tunisien ou au contexte africain.
Packages utilisés
readxl: importer des fichiers Excel.plotly: créer des graphiques interactifs.sf: importer, manipuler et exporter des données géographiques vectorielles.terra: importer, manipuler et exporter des fichiers raster.leaflet: réaliser une carte interactive.geodata: accéder à des jeux de données géographiques de référence dans le monde.rnaturalearth: accéder à des fonds de carte du Monde.wbstats: utiliser les jeux de données de la Banque Mondiale via son API.osmextract: télécharger des données OpenStreetMap.osrm: calculer des temps de parcours routiers via l’engin de routage OSRM.
Si vous n’avez pas ces packages, vous pouvez les installer en exécutant la ligne suivante dans la console.
From scratch
En partant de rien et avec quelques lignes de code on peut créer un jeu de données spatial pour en faire, par exemple, une carte interactive.
Trois vecteurs sont créés, puis transformés en data.frame et enfin convertis en objet spatial avec la fonction st_as_sf() du package sf. Cela permet de créer rapidement une carte interactive avec le package leaflet.
# 3 vecteurs avec nom, longitude et latitude
name <- c("Sentido Bellevue Park", "Université de Sfax")
long <- c(10.579, 10.742)
lat <- c(35.913, 34.736)
# Transformer en data.frame
loc <- data.frame(name, long, lat)
# Transformer en objet spatial
library(sf)
loc <- st_as_sf(loc,
coords = c("long", "lat"),
crs = 4326)
# Visualiser simplement avec leaflet
library(leaflet)
leaflet(loc) |>
addTiles() |>
addCircleMarkers(popup = loc$name) # HôpitauxTableaux de données
Fichiers .csv
Import des données des délégations au format csv. On lui donne le nom de del_df.
del_df <- read.csv("data/tun/don_del.csv", sep = ";", dec = ",")Avec R, il est très important de maîtriser la nature des objets importés ou transformés, et les convertir le cas échéant. La fonction class() nous permet de constater qu’il s’agit d’un data.frame. La fonction str() permet de détailler son contenu.
class(del_df)[1] "data.frame"str(del_df)'data.frame': 264 obs. of 30 variables:
$ del_code : chr "TN.AN.AR" "TN.AN.ET" "TN.AN.KA" "TN.AN.LS" ...
$ del_nom_fr : chr "Ariana Ville" "Cité Ettathamen" "kalaât El Andalous" "Soukra" ...
$ del_nom_ar : chr "أريانة المدينة" "حي التضامن" "قلعة الاندلس" "سكرة" ...
$ gou_code : chr "AN" "AN" "AN" "AN" ...
$ gou_nom : chr "Ariana" "Ariana" "Ariana" "Ariana" ...
$ gou_cap : int 1 0 0 0 0 0 0 1 0 0 ...
$ gou_cap_dist: num 1 8.8 19.8 6.7 8 8.8 15 1 9.5 7.4 ...
$ reg_code : chr "NE" "NE" "NE" "NE" ...
$ reg_nom : chr "Nord-est" "Nord-est" "Nord-est" "Nord-est" ...
$ popto_2004 : int 97687 78311 23045 89151 53752 60896 19404 32329 27977 31792 ...
$ popco_2004 : int 97687 78311 15313 89151 40176 53911 8909 32329 27977 31792 ...
$ immig_2004 : int 16961 5651 728 19129 7258 14053 1298 3963 5831 5508 ...
$ emigr_2004 : int 15426 5245 528 2832 985 1443 723 9381 1135 3991 ...
$ mobil_2004 : int 32387 10896 1256 21961 8243 15496 2021 13344 6966 9499 ...
$ menag_2004 : int 27468 11950 4709 21590 17119 14276 4215 7941 6498 8462 ...
$ ordin_2004 : int 9751 430 188 2979 685 2327 181 834 1176 1963 ...
$ porta_2004 : int 22524 4505 1865 13321 7087 9022 1775 4908 4178 6397 ...
$ telef_2004 : int 18596 3824 829 9262 7019 6339 1176 4002 2755 4747 ...
$ popto_2014 : int 114486 84312 26796 129693 89884 106414 24503 31128 40101 34962 ...
$ popco_2014 : int 114486 84312 18211 129693 58641 94961 11351 31128 40101 34962 ...
$ immig_2014 : int 15637 5028 1104 20786 14400 20128 1997 2392 5855 3957 ...
$ emigr_2014 : int 20448 6752 752 5528 1828 2521 822 8495 1883 3803 ...
$ mobil2014 : int 36085 11780 1856 26314 16228 22649 2819 10887 7738 7760 ...
$ menag_2014 : int 32498 22087 6554 33981 22781 27574 5922 8506 10066 9996 ...
$ ordin_2014 : int 25474 6836 1701 18191 8326 14284 1723 4339 5494 6471 ...
$ porta_2014 : int 32308 21715 6305 33683 22531 27338 5811 8363 9970 9901 ...
$ telef_2014 : int 19942 3496 428 9831 4249 7505 637 2918 3110 4901 ...
$ popto_2010 : num 109500 82922 24367 107829 69247 ...
$ surfa_2010 : num 18.61 3.38 188.21 27.89 22.91 ...
$ idr_2011 : num 0.638 0.386 0.383 0.557 0.466 0.531 0.358 0.489 0.51 0.523 ...On visualise les premières lignes avec la fonction head().
head(del_df, 5) del_code del_nom_fr del_nom_ar gou_code gou_nom gou_cap
1 TN.AN.AR Ariana Ville أريانة المدينة AN Ariana 1
2 TN.AN.ET Cité Ettathamen حي التضامن AN Ariana 0
3 TN.AN.KA kalaât El Andalous قلعة الاندلس AN Ariana 0
4 TN.AN.LS Soukra سكرة AN Ariana 0
5 TN.AN.MN Mnihla المنيهلة AN Ariana 0
gou_cap_dist reg_code reg_nom popto_2004 popco_2004 immig_2004 emigr_2004
1 1.0 NE Nord-est 97687 97687 16961 15426
2 8.8 NE Nord-est 78311 78311 5651 5245
3 19.8 NE Nord-est 23045 15313 728 528
4 6.7 NE Nord-est 89151 89151 19129 2832
5 8.0 NE Nord-est 53752 40176 7258 985
mobil_2004 menag_2004 ordin_2004 porta_2004 telef_2004 popto_2014 popco_2014
1 32387 27468 9751 22524 18596 114486 114486
2 10896 11950 430 4505 3824 84312 84312
3 1256 4709 188 1865 829 26796 18211
4 21961 21590 2979 13321 9262 129693 129693
5 8243 17119 685 7087 7019 89884 58641
immig_2014 emigr_2014 mobil2014 menag_2014 ordin_2014 porta_2014 telef_2014
1 15637 20448 36085 32498 25474 32308 19942
2 5028 6752 11780 22087 6836 21715 3496
3 1104 752 1856 6554 1701 6305 428
4 20786 5528 26314 33981 18191 33683 9831
5 14400 1828 16228 22781 8326 22531 4249
popto_2010 surfa_2010 idr_2011
1 109500 18.612 0.638
2 82922 3.376 0.386
3 24367 188.206 0.383
4 107829 27.895 0.557
5 69247 22.907 0.466Fichiers Excel
La package readxl permet l’import de fichiers Excel. La fonction read_excel a plusieurs arguments utiles pour spécifier en entrée le format des colonnes (texte, numérique), ne pas considérer les n premières lignes du fichier, etc.
Nous importons ici un fichier Excel dérivé du World Population Prospects des Nations Unies, qui met à disposition toute sorte d’indicateurs démographiques (population, natalité, mortalité). On dispose ici des dénombrements de population par tranche d’âge quinquennale de 1950 à 2021 (feuille estimates).
library(readxl)
world <- read_excel("data/world/WPP2022_POP_F02_1_POPULATION_5-YEAR_AGE_GROUPS_BOTH_SEXES.xlsx",
sheet = "Estimates",
skip = 16,
col_types = c(rep("text", 10), rep("numeric", 22))) Manipulation de data.frame
Dans les faits ces fichiers nécessitent généralement d’être reformatés et réorganisés pour pouvoir être interprétés dans un logiciel quel qu’il soit… Commençons par comprendre la structure des données : noms de colonnes, années disponibles etc.
# Conversion en data.frame
world <- data.frame(world)
# Nom des colonnes
colnames(world) [1] "Index"
[2] "Variant"
[3] "Region..subregion..country.or.area.."
[4] "Notes"
[5] "Location.code"
[6] "ISO3.Alpha.code"
[7] "ISO2.Alpha.code"
[8] "SDMX.code.."
[9] "Type"
[10] "Parent.code"
[11] "Year"
[12] "X0.4"
[13] "X5.9"
[14] "X10.14"
[15] "X15.19"
[16] "X20.24"
[17] "X25.29"
[18] "X30.34"
[19] "X35.39"
[20] "X40.44"
[21] "X45.49"
[22] "X50.54"
[23] "X55.59"
[24] "X60.64"
[25] "X65.69"
[26] "X70.74"
[27] "X75.79"
[28] "X80.84"
[29] "X85.89"
[30] "X90.94"
[31] "X95.99"
[32] "X100." # Années disponibles
unique(world$Year) [1] 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961 1962 1963 1964
[16] 1965 1966 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979
[31] 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994
[46] 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009
[61] 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 NAOn peut alors filtrer les dimensions du tableau, calculer des indicateurs, etc. qui pourront être utilisés pour des statistiques de base, des graphiques voire de futures représentations cartographiques.
# Somme de plusieurs colonnes
world$POP_TOT <- rowSums(world[,c(12:32)]) # Population totale
# Renommer une colonne
names(world)[3] <- "Name"
# Créer des indicateurs
# Part des jeunes
world$YOUNG_RT <- (world$X0.4 + world$X5.9 + world$X10.14) / world$POP_TOT * 100
# Part des personnes âgées
world$OLD_RT <- (world$X65.69 + world$X70.74 + world$X75.79 + world$X80.84 +
world$X85.89 + world$X90.94 + world$X95.99) / world$POP_TOT * 100
# Sélection d'une ligne
df <- world[world$Year == 2021,] # Année de référence
# Sélectionner plusieurs lignes
afr <- c("910", "911", "912", "913", "914")
df <- df[df$Parent.code %in% afr,] # Pays africains
# Ordonner selon les valeurs d'une colonne
df <- df[order(df$YOUNG_RT, decreasing = TRUE),]Il existe de nombreuses fonctions pour l’analyse statistique avec R. La plus basique étant probablement summary().
# Résumé stat
summary(df) Index Variant Name Notes
Length:58 Length:58 Length:58 Length:58
Class :character Class :character Class :character Class :character
Mode :character Mode :character Mode :character Mode :character
Location.code ISO3.Alpha.code ISO2.Alpha.code SDMX.code..
Length:58 Length:58 Length:58 Length:58
Class :character Class :character Class :character Class :character
Mode :character Mode :character Mode :character Mode :character
Type Parent.code Year X0.4
Length:58 Length:58 Min. :2021 Min. : 0.22
Class :character Class :character 1st Qu.:2021 1st Qu.: 299.04
Mode :character Mode :character Median :2021 Median : 1971.23
Mean :2021 Mean : 3576.71
3rd Qu.:2021 3rd Qu.: 4280.10
Max. :2021 Max. :34830.80
X5.9 X10.14 X15.19 X20.24
Min. : 0.275 Min. : 0.26 Min. : 0.197 Min. : 0.166
1st Qu.: 287.524 1st Qu.: 260.32 1st Qu.: 231.017 1st Qu.: 218.739
Median : 1799.815 Median : 1596.25 Median : 1340.744 Median : 1073.518
Mean : 3237.939 Mean : 2880.72 Mean : 2481.160 Mean : 2123.874
3rd Qu.: 3812.367 3rd Qu.: 3456.02 3rd Qu.: 3054.270 3rd Qu.: 2727.918
Max. :30567.689 Max. :26974.49 Max. :22929.125 Max. :18806.200
X25.29 X30.34 X35.39 X40.44
Min. : 0.252 Min. : 0.244 Min. : 0.287 Min. : 0.267
1st Qu.: 201.104 1st Qu.: 188.732 1st Qu.: 152.122 1st Qu.: 122.764
Median : 960.347 Median : 880.112 Median : 738.572 Median : 573.226
Mean : 1858.879 Mean : 1643.911 Mean : 1421.126 Mean :1171.518
3rd Qu.: 2325.300 3rd Qu.: 1923.398 3rd Qu.: 1581.452 3rd Qu.:1374.119
Max. :15675.790 Max. :13171.290 Max. :11591.915 Max. :9893.223
X45.49 X50.54 X55.59 X60.64
Min. : 0.358 Min. : 0.486 Min. : 0.452 Min. : 0.496
1st Qu.: 99.267 1st Qu.: 84.229 1st Qu.: 67.378 1st Qu.: 54.439
Median : 469.567 Median : 374.635 Median : 291.823 Median : 233.514
Mean : 938.585 Mean : 770.335 Mean : 616.303 Mean : 473.997
3rd Qu.:1124.716 3rd Qu.: 851.100 3rd Qu.: 667.341 3rd Qu.: 534.762
Max. :7846.818 Max. :6096.186 Max. :4876.240 Max. :3778.750
X65.69 X70.74 X75.79 X80.84
Min. : 0.399 Min. : 0.4345 Min. : 0.3185 Min. : 0.167
1st Qu.: 40.704 1st Qu.: 26.6259 1st Qu.: 16.2689 1st Qu.: 9.041
Median : 156.874 Median : 94.0440 Median : 58.9380 Median : 29.567
Mean : 349.206 Mean : 235.3130 Mean : 140.6908 Mean : 70.866
3rd Qu.: 377.993 3rd Qu.: 237.7533 3rd Qu.: 132.5434 3rd Qu.: 68.634
Max. :2757.823 Max. :1819.8435 Max. :1102.5085 Max. :494.341
X85.89 X90.94 X95.99 X100.
Min. : 0.079 Min. : 0.038 Min. : 0.0075 Min. :0.000000
1st Qu.: 4.053 1st Qu.: 0.968 1st Qu.: 0.0685 1st Qu.:0.004375
Median : 9.991 Median : 2.052 Median : 0.3372 Median :0.031750
Mean : 28.099 Mean : 7.973 Mean : 1.4798 Mean :0.219655
3rd Qu.: 27.383 3rd Qu.: 7.602 3rd Qu.: 1.2139 3rd Qu.:0.152125
Max. :178.652 Max. :81.603 Max. :22.9160 Max. :3.797500
POP_TOT YOUNG_RT OLD_RT
Min. : 5.4 Min. :13.97 Min. : 1.677
1st Qu.: 2388.4 1st Qu.:34.33 1st Qu.: 2.696
Median : 12774.0 Median :40.41 Median : 3.152
Mean : 24028.9 Mean :38.12 Mean : 4.298
3rd Qu.: 28556.3 3rd Qu.:43.45 3rd Qu.: 4.161
Max. :213401.3 Max. :48.90 Max. :26.714 Représentations graphiques
Les possibilités offertes en terme de représentations graphiques sont nombreuses avec R ! Car c’est aussi un langage qui a été créé pour cela.
De base
En une ligne de code on peut créer des représentations graphiques variées, adaptées à la nature statistique des variables à explorer.
par(mar = c(2,2,2,2), mfrow = c(2, 3))
barplot(df$YOUNG_RT, main = "Diagramme en barre")
boxplot(df$YOUNG_RT, main = "Boîtes à moustache")
hist(df$YOUNG_RT, main = "Histogrammes")
hist(df$YOUNG_RT, freq = FALSE, main = "Histogrammes et densité")
lines(density(df$YOUNG_RT), col = "blue")
stripchart(df$YOUNG_RT, method = "jitter", pch = 16,
main = "Diagrammes de dispersion")
plot(data = df, YOUNG_RT ~ OLD_RT, main = "Nuages de points")
Ces graphiques sont paramétrables avec une série d’arguments graphiques.
par(mar = c(4,4,0,4))
# Line plot
tun <- world[world$Name == "Tunisia",]
alg <- world[world$Name == "Algeria",]
mar <- world[world$Name == "Morocco",]
plot(alg$Year, # Abscisses
alg$YOUNG_RT, # Ordonnées
type = "l", # Type lignes
ylim = c(20, 50), # Bornes min/max des ordonnées
cex = .6, # Taille des points
col = "blue", # Couleur de la ligne
cex.lab = 0.7, # Réduit les labels d'un facteur de 0.7
cex.axis = 0.6, # Réduit les labels des graduations d'un facteur de 0.6
xlab = "Années", # Label abscisses
ylab = "Part des jeunes (%) dans la population totale") # Label ordonnées
lines(mar$Year,
mar$YOUNG_RT, # Rajouter une ligne (Maroc)
col = "darkgreen",
type = "l",
cex = .6)
lines(tun$Year, # Et Tunisie
tun$YOUNG_RT,
type = "l",
cex = .6,
col = "red")
# Organisation de la légende
legend(x = "topright",
legend = c("Algérie", "Maroc", "Tunisie"),
col = c("blue", "darkgreen", "red"), lty = 1,
cex = .8)
Graphiques interactifs
Le package plotly (Sievert et al. 2024) permet d’intégrer une dimension interactive aux représentations graphiques. Pour utiliser les fonctions de ce package, il faut bien avoir en tête le format de données attendues (format long).
library(plotly)
# Importer jeu de données d'exemple de gapminder
df <- read.csv("data/world/gapminder.csv")
head(df) country continent year lifeExp pop gdpPercap
1 Afghanistan Asia 1952 28.801 8425333 779.4453
2 Afghanistan Asia 1957 30.332 9240934 820.8530
3 Afghanistan Asia 1962 31.997 10267083 853.1007
4 Afghanistan Asia 1967 34.020 11537966 836.1971
5 Afghanistan Asia 1972 36.088 13079460 739.9811
6 Afghanistan Asia 1977 38.438 14880372 786.1134# Créer un graphique de type "Multiple Trace Animations"
fig <- df %>%
plot_ly(
x = ~gdpPercap,
y = ~lifeExp,
size = ~pop,
color = ~continent,
frame = ~year,
text = ~country,
hoverinfo = "text",
type = 'scatter',
mode = 'markers'
)
# Paramétrer l'échelle, les labels, etc.
fig <- fig %>% layout(
xaxis = list(type = "log", title = "PIB par habitant"),
yaxis = list(title = "Espérance de vie")
)
fig
Pour aller plus loin
Consulter r-graph-gallery qui présente les possibilités graphiques les plus courantes, en langage de base R. L’exploration de la syntaxe ggplotet son package éponyme ggplot2 peut être utile pour la recherche de représentations graphiques avancées. Cette syntaxe est un peu différente du langage de base R et repose sur les principes de la “grammaire des graphiques”. Plusieurs manuels très bien construits permettent de rentrer dans l’univers de la visualisation de données avec ggplot2, comme Modern Data Visualization with R (Kabacoff 2023)
L’intérêt de la programmation pour manipuler des données
Les fournisseurs de données institutionnels distribuent parfois des tableaux de données peu adaptés à l’intégration dans un logiciel quel qu’il soit. La mise en forme des données dans un langage de programmation permet :
- d’éviter d’avoir à manipuler ces fichiers manuellement, source d’erreurs.
- d’avoir un façon générique et rapide d’importer un grand nombre de tableaux de données (et éviter de perdre de nombreuses heures de travail).
Voici un exemple adapté aux données de chômage par gouvernorats tunisiens produit par l’Institut National de la Statistique tunisien (INS). Le fichier brut ressemble à ceci :
Un court bloc de code permet de réorganiser le fichier comme désiré, et pourrait être utilement étendu aux autres fichiers distribués par l’INS.
# Import
df <- read_excel("data/tun/chomage.xls", # Chemin du fichier
sheet = "Sheet1", # Nom de la feuille Excel
skip = 3, # Retirer les trois premières lignes
col_types = c(rep("text", 3), rep("numeric", 6))) # Format des colonnes
# Table de passage Recensement / codes internationaux
gouv <- read.csv("data/tun/table_GADM_TUN.csv")
# Reformatage
df <- data.frame(df)
df[,1:2] <- NULL # Retirer les deux premières colonnes
# Noms de colonnes
cols <- c("CHOM_INS_NON_", "CHOM_INS_PRI_", "CHOM_INS_SEC_", "CHOM_INS_SUP_",
"CHOM_INS_UNK_", "CHOM_INS_TOT_")
colnames(df)[1] <- "id_TUN"
# Chômeurs
# Hommes
df1 <- df[c(1:25),]
names(df1)[2:length(df1)] <- paste0(cols, "M")
gouv <- merge(gouv, df1, by = "id_TUN", all.x = TRUE)
# Femmes
df2 <- df[c(26:50),]
names(df2)[2:length(df2)] <- paste0(cols, "F")
gouv <- merge(gouv, df2, by = "id_TUN", all.x = TRUE)
# Total
df3 <- df[c(51:75),]
names(df3)[2:length(df3)] <- paste0(cols, "T")
gouv <- merge(gouv, df3, by = "id_TUN", all.x = TRUE)
head(gouv) id_TUN id_GADM CHOM_INS_NON_M CHOM_INS_PRI_M CHOM_INS_SEC_M CHOM_INS_SUP_M
1 ARIANA TN.AN 484 2711 6876 2979
2 BEJA TN.BJ 1211 2701 4577 1399
3 BEN AROUS TN.BA 400 3496 8869 3512
4 BIZERTE TN.BZ 1276 4592 8087 2568
5 GABES TN.GB 305 2465 4890 2620
6 GAFSA TN.GF 488 3031 7571 3576
CHOM_INS_UNK_M CHOM_INS_TOT_M CHOM_INS_NON_F CHOM_INS_PRI_F CHOM_INS_SEC_F
1 69 13119 469 1459 4674
2 77 9965 1025 1319 2728
3 81 16358 696 2024 5960
4 74 16597 698 1775 4063
5 100 10380 340 1155 3007
6 67 14733 531 1711 4516
CHOM_INS_SUP_F CHOM_INS_UNK_F CHOM_INS_TOT_F CHOM_INS_NON_T CHOM_INS_PRI_T
1 6610 44 13256 953 4170
2 2718 51 7841 2236 4020
3 6977 37 15694 1096 5520
4 5311 64 11911 1974 6367
5 7235 100 11837 645 3620
6 7639 80 14477 1019 4742
CHOM_INS_SEC_T CHOM_INS_SUP_T CHOM_INS_UNK_T CHOM_INS_TOT_T
1 11550 9589 113 26375
2 7305 4117 128 17806
3 14829 10489 118 32052
4 12150 7879 138 28508
5 7897 9855 200 22217
6 12087 11215 147 29210Une fois mis en forme, on peut réaliser un graphique sans avoir modifié le fichier initial provenant de l’INS.
# 2 graphiques par ligne
par(mfrow = c(1,2), mar = c(2,2,2,2))
# Boxplot hommes
boxplot(gouv$CHOM_INS_NON_M, # Chômeurs hommes par niveau d'instruction
gouv$CHOM_INS_PRI_M,
gouv$CHOM_INS_SEC_M,
gouv$CHOM_INS_SUP_M,
ylim = c(0, 12000), # Bornes min/max de l'axe des ordonnées
main = "Hommes", # Titre plot
names = c("Rien", "Primaire", "Secondaire", "Tertiaire"), # Labels (X)
ylab = "Nombre de chômeurs, 2014", # Label (Y)
col = "#adcaf7", # Couleur des box-plots
cex.axis = .6, # Taille des labels des axes (réduit de 70 %)
cex.title = .6) # Taille du label du titre (réduit de 70 %)
# Boxplot femmes
boxplot(gouv$CHOM_INS_NON_F,
gouv$CHOM_INS_PRI_F,
gouv$CHOM_INS_SEC_F,
gouv$CHOM_INS_SUP_F,
ylim = c(0, 12000),
main = "Femmes",
names = c("Rien", "Primaire", "Secondaire", "Tertiaire"),
col = "#ed9fb0",
cex.axis = .6,
cex.title = .6)
Export
write.csv() exporte un data.frame selon un chemin spécifié.
write.csv(x = gouv, # Objet à exporter
dsn = "data/tun/gouv_chom.csv", # Chemin de fichier
row.names = FALSE) # Pour retirer les numéros de ligneDonnées vectorielles
Le package sf permet d’importer et manipuler des couches géographiques vectorielles (points, lignes, polygones).
Import
Le geopackage est un format de données ouvert qui permet de stocker plusieurs couches géographiques dans un même fichier. La fonction st_layers() permet d’avoir un aperçu des couches présentes dans un fichier geopackage.
library(sf)
st_layers("data/tun/geom/tun_admin.gpkg")Driver: GPKG
Available layers:
layer_name geometry_type features fields crs_name
1 delegation Multi Polygon 266 11 ETRS89-extended / LAEA Europe
2 gouvernorat Multi Polygon 24 4 ETRS89-extended / LAEA Europe
3 region Multi Polygon 6 2 ETRS89-extended / LAEA EuropeImporter des couches géographiques avec st_read(). Il s’agit ici des couches géographiques créées il y a quelques années par l’UAR RIATE.
del <- st_read("data/tun/geom/tun_admin.gpkg", layer = "delegation", quiet = TRUE)
gouv <- st_read("data/tun/geom/tun_admin.gpkg", layer = "gouvernorat", quiet = TRUE)
reg <- st_read("data/tun/geom/tun_admin.gpkg", layer = "region", quiet = TRUE)Il s’agit bien d’objets sf.
class(del)[1] "sf" "data.frame"La fonction st_read() peut aussi être employée pour des formats de fichiers .geojson, .shapefiles, etc.
reg <- st_read("data/tun/geom/map_reg.geojson", quiet = TRUE) Les objets sf sont des data.frame dont l’une des colonnes contient des géométries. Cette colonne est de la classe sfc (simple feature column) et chaque individu de la colonne est un sfg (simple feature geometry).
Ce format est très pratique dans la mesure où les données et les géométries sont intrinsèquement liées dans un même objet. Des informations additionnelles, propres aux couches d’informations géographiques sont aussi restituées (type de géométrie, emprise spatiale, système de coordonnées de référence).
head(del)Simple feature collection with 6 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 4184398 ymin: 1282883 xmax: 4380989 ymax: 1539573
Projected CRS: ETRS89-extended / LAEA Europe
del_code del_nom_fr del_code_riate del_code_ins del_nom_ar gou_code
1 TN.SF.AG Agareb TS2146 TN347 عقارب TN.SF
2 TN.JE.AD Ain Drahem TS1224 TN225 عين دراهم TN.JE
3 TN.SS.AK Akouda TS2115 TN316 اكودة TN.SS
4 TN.KR.AL Alaa TS2216 TN417 العلاء TN.KR
5 TN.BJ.AM Amdoun TS1212 TN213 عمدون TN.BJ
6 TN.AN.AR Ariana Ville TS1120 TN121 أريانة المدينة TN.AN
gou_nom gou_cap gou_cap_dist reg_code reg_nom
1 Sfax 0 30.1 CE Centre-est
2 Jendouba 0 38.5 NO Nord-ouest
3 Sousse 0 14.3 CE Centre-est
4 Kairouan 0 48.6 CO Centre-ouest
5 Beja 0 18.9 NO Nord-ouest
6 Ariana 1 1.0 NE Nord-est
geom
1 MULTIPOLYGON (((4375308 130...
2 MULTIPOLYGON (((4224769 152...
3 MULTIPOLYGON (((4373533 142...
4 MULTIPOLYGON (((4306097 141...
5 MULTIPOLYGON (((4243838 153...
6 MULTIPOLYGON (((4338370 153...Aperçu des variables avec plot() :
plot(del)
Afficher juste les géométries:
# 3 cartes par ligne
par(mfrow = c(1,3), mar = c(2,2,2,2))
# Délégations
plot(del$geom, # Géométries uniquement
col = "peachpuff", # Couleur du fond
border = "white", # Couleur de bordure
main = "Délégations") # Titre
# Gouvernorats
plot(gouv$geom,
col = "peachpuff",
border = "white",
main = "Gouvernorats")
# "Régions"
plot(reg$geom,
col = "peachpuff",
border = "white",
main = "Régions")
Jointures attributaires
Nous pouvons joindre un data.frame à un objet sf en utilisant la fonction merge() et en s’appuyant sur des identifiants communs aux deux objets.
Attention à l’ordre des arguments, l’objet retourné sera du même type que x. Il n’est pas possible de faire une jointure attributaire en utilisant deux objets sf.
del <- merge(
x = del[,"del_code"], # L'objet sf (seulement le champ del_code)
y = del_df, # le data.frame
by.x = "del_code", # identifiant dans x
by.y = "del_code", # identifiant dans y
all.x = TRUE # conserver toutes les lignes
)
head(del)Simple feature collection with 6 features and 30 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 4321914 ymin: 1528813 xmax: 4346090 ymax: 1562156
Projected CRS: ETRS89-extended / LAEA Europe
del_code del_nom_fr del_nom_ar gou_code gou_nom gou_cap
1 TN.AN.AR Ariana Ville أريانة المدينة AN Ariana 1
2 TN.AN.ET Cité Ettathamen حي التضامن AN Ariana 0
3 TN.AN.KA kalaât El Andalous قلعة الاندلس AN Ariana 0
4 TN.AN.LS Soukra سكرة AN Ariana 0
5 TN.AN.MN Mnihla المنيهلة AN Ariana 0
6 TN.AN.RA Raoued رواد AN Ariana 0
gou_cap_dist reg_code reg_nom popto_2004 popco_2004 immig_2004 emigr_2004
1 1.0 NE Nord-est 97687 97687 16961 15426
2 8.8 NE Nord-est 78311 78311 5651 5245
3 19.8 NE Nord-est 23045 15313 728 528
4 6.7 NE Nord-est 89151 89151 19129 2832
5 8.0 NE Nord-est 53752 40176 7258 985
6 8.8 NE Nord-est 60896 53911 14053 1443
mobil_2004 menag_2004 ordin_2004 porta_2004 telef_2004 popto_2014 popco_2014
1 32387 27468 9751 22524 18596 114486 114486
2 10896 11950 430 4505 3824 84312 84312
3 1256 4709 188 1865 829 26796 18211
4 21961 21590 2979 13321 9262 129693 129693
5 8243 17119 685 7087 7019 89884 58641
6 15496 14276 2327 9022 6339 106414 94961
immig_2014 emigr_2014 mobil2014 menag_2014 ordin_2014 porta_2014 telef_2014
1 15637 20448 36085 32498 25474 32308 19942
2 5028 6752 11780 22087 6836 21715 3496
3 1104 752 1856 6554 1701 6305 428
4 20786 5528 26314 33981 18191 33683 9831
5 14400 1828 16228 22781 8326 22531 4249
6 20128 2521 22649 27574 14284 27338 7505
popto_2010 surfa_2010 idr_2011 geometry
1 109500 18.612 0.638 MULTIPOLYGON (((4338370 153...
2 82922 3.376 0.386 MULTIPOLYGON (((4330561 153...
3 24367 188.206 0.383 MULTIPOLYGON (((4333853 156...
4 107829 27.895 0.557 MULTIPOLYGON (((4342368 153...
5 69247 22.907 0.466 MULTIPOLYGON (((4332840 153...
6 77419 57.080 0.531 MULTIPOLYGON (((4338691 154...plot(del[,"idr_2011"])
Sélectionner des lignes, des colonnes
Les objets sf sont des data.frame, on peut donc sélectionner leur lignes et leur colonnes de la même manière.
# Sélection de lignes
sou <- del[del$gou_nom == "Sousse",]
# Sélection de colonnes
sou <- sou[,"idr_2011"]
# Ne conserver que les lignes avec une valeur
sou <- sou[!is.na(sou$idr_2011),]
plot(sou[,"idr_2011"])
Export d’une couche géographique
st_write("data/sousse_deleg.geojson")Données raster
Le package terra (Hijmans 2024) permet d’importer et d’exporter des fichiers raster. Il repose sur la bibliothèque GDAL (GDAL/OGR contributors 2022) qui permet de lire et de traiter un très grand nombre de format d’images géographiques.
Import
La fonction rast() permet de créer et/ou d’importer des données raster. Nous importons ici un jeu de données au format .tif créé et distribué par WorldPop qui porte sur une estimation de la population tunisienne dans une résolution de 100 mètres.
library(terra)
pop <- rast("data/tun/raster/tun_ppp_2020_UNadj.tif")Ce sont des objets de type SpatRaster.
popclass : SpatRaster
dimensions : 8788, 4881, 1 (nrow, ncol, nlyr)
resolution : 0.0008333333, 0.0008333333 (x, y)
extent : 7.530417, 11.59792, 30.23708, 37.56042 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : tun_ppp_2020_UNadj.tif
name : tun_ppp_2020_UNadj La fonction summary() est toujours utile pour un résumé statistique des cellules. Vu le nombre important de cellules, ce résumé est effectué sur un échantillon.
summary(pop) tun_ppp_2020_UNadj
Min. : 0.00
1st Qu.: 0.01
Median : 0.09
Mean : 0.55
3rd Qu.: 0.25
Max. :234.75
NA's :49414 Affichage
Par défaut et comme montré précédemment, la fonction plot() renvoie une légende continue sur les représentations cartographiques. On peut choisir de discrétiser cette information pour représenter des classes de valeur associée à une palette de couleurs personnalisée.
par(mfrow = c(1,2))
# Graphique en échelle continue
plot(pop)
# Avec discrétisation et paramétrage des couleurs
cuts <- c(0, 1, 2, 4, 8, 16, 32, 64, 235)
cols <- colorRampPalette(c("yellow", "darkgoldenrod1", "brown1"))
plot(pop, breaks = cuts, col = cols(8))
Modification de la zone d’étude
Le découpage d’un raster en fonction de l’étendue d’un autre objet, SpatVector ou SpatRaster, est réalisable avec la fonction crop(). Les deux couches de données doivent être dans la même projection.
# Transformer en WGS 84
sou <- st_transform(sou, 4326)
# Crop avec la délégation de Sousse
pop_sou <- crop(pop, sou)
# Plot
plot(pop_sou, breaks = cuts, col = cols(8))
plot(sou$geometry, col = NA, add = TRUE)
Export
La fonction writeRaster() permet d’enregistrer un objet SpatRaster.
writeRaster(x = pop_sou, filename = "data/tun/raster/pop_sou.tif")Packages de données spatiales
De nombreux packages distribuent des données géographiques. Ils interfacent souvent des API qui permettent d’interroger des données mises à disposition sur le Web, directement avec R.
L’intérêt des API est de se connecter directement aux bases de données des fournisseurs de données, et de disposer des dernières mises à jour des données.
Principaux packages
Au niveau international
rnaturalearth(Massicotte et South 2023) : permet de récupérer les données cartographiques Natural Earth.
cshapes(Weidmann, Schvitz, et Girardin 2021) : met à disposition les frontières nationales, de 1886 à aujourd’hui.
geonames(Rowlingson 2019) : permet d’interroger la BD geonames, qui fournit notamment des localisations.wbstats(Piburn 2020) etWDI(Arel-Bundock 2022) : donnent accès aux données et statistiques de la Banque mondiale.
sen2r(Ranghetti et al. 2020) : permet de télécharger et prétraiter automatiquement les données du satellite Sentinel-2.
MODIStsp(Busetto et Ranghetti 2016) : permet de trouver, télécharger et traiter des images MODIS.geodata(Hijmans et al. 2023) : fournit un accès à des données sur le climat, l’altitude, le sol, la présence d’espèces et les limites administratives.elevatr(Hollister et al. 2023) : donne accès à des données d’élévation mises à disposition par Amazon Web Services Terrain Tiles, l’Open Topography Global Datasets API et l’USGS Elevation Point Query Service.rgee(Aybar 2023) : permet d’utiliser l’API de Google Earth Engine, catalogue de données publiques et infrastructure de calcul pour les images satellites.
nasapower(Sparks 2018) : API client NASA (prévision des ressources énergétiques mondiales, météorologie, énergie solaire de surface et climatologie).geoknife(Read et al. 2015) : permet de traiter (en ligne) des données matricielles volumineuses issues du Geo Data Portal de l’U.S. Geological Survey.rdhs(Watson, FitzJohn, et Eaton 2019) : API client et gestions de données de l’enquête démographique et de santé (DHS).giscoR(Hernangómez 2024a) : permet de télécharger des données cartographiques mondiales et européennes de la BD GISCO d’Eurostat (système d’information géographique de la Commission).eurostat(Lahti et al. 2017) : permet de télécharger des données de la BD Eurostat.
Au niveau national
- Brésil
geobr(Pereira et Goncalves 2024) : fournit un accès facile aux séries de données spatiales officielles du Brésil pour différentes années et découpage administratifs.
- Chili -
chilemapas(Vargas 2022) : donne accès aux divisions politiques et administratives du Chili. - Espagne
mapSpain(Hernangómez 2024b) : propose les limites administratives de l’Espagne à plusieurs niveaux (Communautés autonomes, Provinces, Municipalités), ainsi que des tuiles.
- États-Unis
tidycensus(Walker et Herman 2024) : permet de charger des données et géométries du recensement américain en formatsfettidyversetigris(Walker 2024b) : donne accès aux éléments cartographiques fournis par le US Census Bureau TIGER, y compris les limites cartographiques, les routes et l’eau.FedData(Bocinsky 2024) : automatise le téléchargement de données géospatiales disponibles à partir de plusieurs sources de données fédérées.acs(Glenn 2019) : permet de télécharger et manipuler les données de l’American Community Survey et les données décennales du recensement des États-Unis.censusapi(Recht 2022) : wrapper pour les API du Census Bureau des États-Unis.idbr(Walker 2024a) : interface avec l’API de la base de données internationale du US Census Bureau.
ipumsr(Greg Freedman Ellis, Derek Burk, et Finn Roberts 2024) : Permet d’importer des données de recensement, d’enquête et géographiques fournies par l’IPUMS.totalcensus(Li 2021) : permet d’extraire les données du recensement décennal et de l’American Community Survey au niveaux des block, block group et tract.
- Finland
mapsFinland(Haukka 2023) : donne un accès à des cartes et données concernant la Finland.
- France
happign(Carteron 2023) : accès à certaines données de l’IGN.insee(Leclerc 2022) : pour télécharger facilement les données de la base BDM (Banque de Données Macroéconomiques) de l’INSEE.
- Pologne
rgugik(Dyba et Nowosad 2021) : permet l’acquisition automatique de données ouvertes à partir des ressources du Bureau central polonais de la géodésie et de la cartographie (Główny Urząd Geodezji i Kartografii ).
- Uruguay
geouy(Detomasi 2023) : permet le chargement d’informations géographiques sur l’Uruguay.
- …
geodata
Ce package facilite l’accès à des données géographique de référence sur le climat, la couverture du sol, les limites administratives et plusieurs autres jeux de données de référence au niveau mondial.
Les couches géographiques inclues dans le package
geodata
| Function | Description |
|---|---|
bio_oracle |
Marine data from Bio-Oracle |
cmip6_world, cmip6_tile |
Downscaled and calibrated CMIP6 projected future climate data |
country_codes |
Country codes |
crop_calendar_sacks |
Crop calendar data by Sacks et al |
crop_monfreda |
Crop area and yield data for 175 crops by Monfreda et al. |
crop_spam |
MapSPAM crop data (area, yield, value) |
cropland |
Cropland density for the world derived from different sources (ESA, GLAD, QED) |
elevation_3s, elevation_30s, elevation_global |
Elevation data |
gadm, world |
Administrative boundaries for any country, or the entire world from GADM |
landcover |
Landcover data derived from ESA WorldCover |
footprint |
Human footprint data from the Last of the Wild project |
osm |
OpenStreetMap data by country (places and roads) |
population |
Population density Gridded Population of the World |
soil_af |
Chemical and physical soil properties data for Africa for different soil depths |
soil_af_elements |
Connect to or download chemical soil element concentration (for the 0-30 cm topsoil) data for Africa |
soil_af_water |
Physical soil properties data for Africa for water balance computation |
soil_af_isda |
Soil data for Africa derived from the iDSA data set |
soil_world_vsi |
Virtually connect to the global SoilGrids data |
soil_world |
Global soils data from SoilGrids |
sp_occurrence |
Species occurrence records from the Global Biodiversity Information Facility |
travel_time |
Travel time to and from cities and ports by Nelson et al. |
worldclim_global, worldclim_country, worldclim_tile |
WorldClim glocal climate data |
Extraction des couches géographiques d’altitude et de température en Tunisie.
library(geodata)
elev <- elevation_30s(country = "TUN", path = tempdir())
temp <- worldclim_country(country = "Tunisia",
res = 10,
var = "tavg",
path = tempdir())Ce sont des objets de type SpatRaster.
class(temp)[1] "SpatRaster"
attr(,"package")
[1] "terra"class(elev)[1] "SpatRaster"
attr(,"package")
[1] "terra"# Les afficher
par(mfrow = c(1,2))
# Altitude
cols <- colorRampPalette(c("#31ad37", "#f5f752", "#fca330", "#9c5903"))
plot(elev, main = "Altitude", col = cols(50))
# Température
cols <- colorRampPalette(c("#4575b4", "#91bfdb", "#e0f3f8",
"#fee090", "#fc8d59", "#d73027"))
plot(temp$TUN_wc2.1_30s_tavg_5, main = "Températures, Mai (1970-2000)",
col = cols(50))
Natural Earth et Banque Mondiale.
Objectif
Nous cherchons ici à créer une carte des émissions de CO² en Afrique en utilisant les données de la banque mondiale et le fond de carte Natural Earth. Rien besoin de télécharger, mais il faut une connexion internet.
Plusieurs fonds de carte des pays du Monde à différents niveaux de généralisation cartographique sont disponibles avec la package rnaturalearth. Nous sélectionnons uniquement les pays africains.
library(rnaturalearth)
# Import des pays
country <- ne_countries(type = "countries", # pays
scale = "small", # niveau de généralisation
returnclass = "sf") # objet retourné
# Conversion en projection Mercator
country <- st_transform(country, crs = "EPSG:3857")
# Si on s'intéresse à l'Afrique (modèle carto)
afr <- country[country$continent == "Africa",]Le package wbstatspermet d’interroger l’API de la base de données de la Banque Mondiale. On peut faire une recherche pour trouver le nom des tables qui répondent à une requête par mots-clés avec la fonction wb_search().
library(wbstats)
wb_search("CO2 emissions")# A tibble: 20 × 3
indicator_id indicator indicator_desc
<chr> <chr> <chr>
1 EN.ATM.CO2E.CP.KT CO2 emissions from cement production (th… "Carbon dioxi…
2 EN.ATM.CO2E.FF.KT CO2 emissions from fossil-fuels, total (… "Fossil fuel …
3 EN.ATM.CO2E.FF.ZS CO2 emissions from fossil-fuels (% of to… "Fossil fuel …
4 EN.ATM.CO2E.GF.KT CO2 emissions from gaseous fuel consumpt… "Carbon dioxi…
5 EN.ATM.CO2E.GF.ZS CO2 emissions from gaseous fuel consumpt… "Carbon dioxi…
6 EN.ATM.CO2E.GL.KT CO2 emissions from gas flaring (thousand… "Carbon dioxi…
7 EN.ATM.CO2E.KD.GD CO2 emissions (kg per 2010 US$ of GDP) "Carbon dioxi…
8 EN.ATM.CO2E.KT CO2 emissions (kt) "Carbon dioxi…
9 EN.ATM.CO2E.LF.KT CO2 emissions from liquid fuel consumpti… "Carbon dioxi…
10 EN.ATM.CO2E.LF.ZS CO2 emissions from liquid fuel consumpti… "Carbon dioxi…
11 EN.ATM.CO2E.PC CO2 emissions (metric tons per capita) "Carbon dioxi…
12 EN.ATM.CO2E.PP.GD CO2 emissions (kg per PPP $ of GDP) "Carbon dioxi…
13 EN.ATM.CO2E.PP.GD.KD CO2 emissions (kg per 2017 PPP $ of GDP) "Carbon dioxi…
14 EN.ATM.CO2E.SF.KT CO2 emissions from solid fuel consumptio… "Carbon dioxi…
15 EN.ATM.CO2E.SF.ZS CO2 emissions from solid fuel consumptio… "Carbon dioxi…
16 EN.CO2.BLDG.ZS CO2 emissions from residential buildings… "CO2 emission…
17 EN.CO2.ETOT.ZS CO2 emissions from electricity and heat … "CO2 emission…
18 EN.CO2.MANF.ZS CO2 emissions from manufacturing industr… "CO2 emission…
19 EN.CO2.OTHX.ZS CO2 emissions from other sectors, exclud… "CO2 emission…
20 EN.CO2.TRAN.ZS CO2 emissions from transport (% of total… "CO2 emission…La fonction wb_data() extrait les tables identifiées pour pour construire l’indicateur d’émissions de CO² par habitant. La série statistique est disponible de 1960 à 2023.
Avertissement
Chaque package reposant sur une API a son fonctionnement propre. Il est donc conseillé de consulter attentivement la documentation associée à ces packages (souvent précise) pour apprendre à les utiliser.
# Sélection des indicateurs
my_indicators = c("pop" = "SP.POP.TOTL",
"co2" = "EN.ATM.CO2E.KT")
# Interroger l'API
wb <- wb_data(my_indicators, return_wide = FALSE)
head(wb)# A tibble: 6 × 11
indicator_id indicator iso2c iso3c country date value unit obs_status
<chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr>
1 SP.POP.TOTL Population, t… AF AFG Afghan… 2023 NA <NA> <NA>
2 SP.POP.TOTL Population, t… AF AFG Afghan… 2022 4.11e7 <NA> <NA>
3 SP.POP.TOTL Population, t… AF AFG Afghan… 2021 4.01e7 <NA> <NA>
4 SP.POP.TOTL Population, t… AF AFG Afghan… 2020 3.90e7 <NA> <NA>
5 SP.POP.TOTL Population, t… AF AFG Afghan… 2019 3.78e7 <NA> <NA>
6 SP.POP.TOTL Population, t… AF AFG Afghan… 2018 3.67e7 <NA> <NA>
# ℹ 2 more variables: footnote <chr>, last_updated <date># Années disponibles
unique(wb$date) [1] 2023 2022 2021 2020 2019 2018 2017 2016 2015 2014 2013 2012 2011 2010 2009
[16] 2008 2007 2006 2005 2004 2003 2002 2001 2000 1999 1998 1997 1996 1995 1994
[31] 1993 1992 1991 1990 1989 1988 1987 1986 1985 1984 1983 1982 1981 1980 1979
[46] 1978 1977 1976 1975 1974 1973 1972 1971 1970 1969 1968 1967 1966 1965 1964
[61] 1963 1962 1961 1960Nous gardons la dernière année pour laquelle les émissions de CO² sont disponibles 2020. La fonction reshape() est utilisée pour transformer le data.frame dans un format attendu (long vers large).
# Sélectionner la dernière année disponible et les colonnes utiles
wb <- wb[wb$date %in% 2020,]
wb <- wb[,c("iso3c", "indicator_id", "value")]
wb <- data.frame(wb)
# Formatage (format long + concaténation code indicateur + année ref
wb <- reshape(wb, idvar = "iso3c", timevar = "indicator_id", direction = "wide")
# Renommage des colonnes et création de l'indicateurs émissions par hab
names(wb)[2:3] <- c("co2_2020", "pop_2020")
wb$co2_hab_2020 <- wb$co2_2020 * 1000 / wb$pop_2020La jointure attributaire avec le fond de carte précédemment importé avec rnaturalearth permet ainsi de constater que tout est en place pour procéder à une représentation cartographique du phénomène.
# Jointure attributaire
afr <- merge(afr[,c("adm0_a3", "name_fr")], wb,
by.x = "adm0_a3", by.y = "iso3c", all.x = TRUE)
# Carte
plot(afr[,"co2_hab_2020"])
OpenStreetMap
OpenStreetMap (OSM) est un projet de cartographie participative qui a pour but de constituer une base de données géographiques libre à l’échelle mondiale. OSM vous permet de voir, modifier et utiliser des données géographiques dans le monde entier.
Plusieurs packages permettre d’extraire, interroger et visualiser des données issues d’OSM :
osmextract(Gilardi et Lovelace 2023) : permet d’importer des données OSM.
osmdata(Mark Padgham et al. 2017) : pour télécharger et utiliser les données d’OSM, en utilisant l’API Overpass turbo.maptiles(Giraud 2024) : Ce package télécharge, compose et affiche des tuiles à partir d’un grand nombre de fournisseurs (OpenStreetMap, Stadia, Esri, CARTO ou Thunderforest…).
nominatimlite(Hernangómez et Lacko 2024) ettidygeocoder(Cambon et al. 2021) : Des géocodeurs qui utilisent OSM, notamment.
Objectif
Dans un premier temps nous cherchons à extraire d’OSM les géométries suivantes en Tunisie :
- Découpages administratifs : Délégations et secteurs.
- Équipements : objets décrits avec une clé OSM amenity (équipements utiles et importants).
A partir de ces données extraites, nous proposons ensuite de calculer les temps de parcours routiers des centres de secteur vers l’hôpital ou la clinique la plus accessible.
Extraction et préparation des données
Le package osmextract (Gilardi et Lovelace 2023) permet d’extraire des données depuis une base de données OSM directement et travailler en local sur des volumes de données très importants et ainsi d’éviter de surcharger un serveur Overpass turbo (utilisé par le package osmdata).
Découpages administratifs
La fonction oe_get() permet de télécharger un extrait de la base de données OSM pour une zone particulière et un type d’objet géographique. L’argument place correspond au nom du fichier *.pbf accessible sur le site Geofabrik. L’argument extra_tag permet de sélectionner les objets de la base OSM correspondant à une clé particulière (se référer à la documentation d’OSM pour choisir les clés).
On commence par extraire tous les polygones et les points OSM inclus en Tunisie. L’import d’extra_tags permet d’obtenir des informations supplémentaires qui serviront à filtrer la base de données.
library(osmextract)
osm_poly <- oe_get(place = "Tunisia",
layer = "multipolygons",
extra_tags = c("amenity", "ref:tn:hasc_2", "ref:tn:codegeo",
"name:fr", "name:ar"))
osm_pt <- oe_get(place = "Tunisia",
layer = "points",
extra_tags = c("amenity", "name:fr", "name:ar"))On peut ensuite filtrer l’extraction en fonction des objets / champs d’intérêt.
# Délégations et secteurs
admin <- osm_poly[!is.na(osm_poly$admin_level),] # Retirer pas d'attribut de niveau hiérarchique
deleg <- admin[admin$admin_level == 5,] # Délégations
sect <- admin[admin$admin_level == 6,] # Secteurs
# Ne garder que les champs utiles
deleg <- deleg[,c("osm_id", "ref_tn_codegeo", "ref_tn_hasc_2", "name_ar",
"name_fr", "admin_level")]
sect <- sect[,c("osm_id", "ref_tn_codegeo", "name_ar",
"name_fr", "admin_level")]Ce bloc de code, dont on ne détaillera pas le contenu ici, permet d’harmoniser le nom des champs entre les couches géographiques et de disposer des délégations et gouvernorats d’appartenance des couches géographiques des secteurs et délégations.
# Délégation d'appartenance du secteur
sect_pt <- st_centroid(sect) # Centroid du secteur
sect_pt <- st_intersection(sect_pt, deleg) # Intersection avec couche délégation
sect_pt <- st_set_geometry(sect_pt, NULL) # Retirer géométries
sect <- merge(sect, # Enrichir les secteurs du code d'appartenance de la délégation
sect_pt[,c("osm_id", "ref_tn_hasc_2")],
by = "osm_id",
all.x = TRUE)
sect <- sect[!duplicated(sect$osm_id),] # Retirer valeurs dupliquées
# Renommer colonnes
names(sect)[2] <- "id_tn"
names(sect)[6] <- "id_hasc_deleg"
sect$id_hasc_gouv <- substr(sect$id_hasc_deleg, 1, 5) # Gouvernorat d'appartenance
names(sect)[2,6] <- (c("tn_codegeo", "int_codegeo"))
names(deleg)[2:3] <- c("id_tn", "id_hasc_deleg")
deleg$id_hasc_gouv <- substr(deleg$id_hasc_deleg, 1, 5)Équipements
Les objets décrits par la clé amenity sont de nature hétérogènes (points, polygones). Ce bloc de code extrait les centroides des polygones et les associent à la couche de points initiale.
Ces couches géographiques sont finalement exportées dans un geopackage tun_osm.gpkg.
# Consolidation des géométries
sel_poly <- osm_poly[!is.na(osm_poly$amenity),] # Retirer les objets qui n'ont pas de clé amenity
sel_pt <- osm_pt[!is.na(osm_pt$amenity),]
sel_poly <- st_make_valid(sel_poly) # Consolider les géométries des polygones
sel_pt2 <- st_centroid(sel_poly) # extraire le centroide
cols <- intersect(names(sel_pt), names(sel_pt2)) # Garder les colonnes identiques
sel_pt <- rbind(sel_pt[,cols, sel_pt2[,cols]]) # Combiner points et polygones
sel_pt <- sel_pt[,c("osm_id", "amenity", "name_ar", "name_fr")]
# Exporter les couche ainsi créées
st_write(sel_pt, "data/tun/geom/tun_osm.gpkg", layer = "poi")
st_write(deleg, "data/tun/geom/tun_osm.gpkg", layer = "deleg")
st_write(sect, "data/tun/geom/tun_osm.gpkg", layer = "sect")
Pour aller plus loin
Ce document (Giraud 2017) montre comment créer un fond de carte avec des données OSM, extraire des objets d’intérêt (bars et restaurants) avec le package osmdata. Il propose également plusieurs pistes cartographiques pour manipuler et visualiser ces données. Notons simplement que la cartographie utilise les fonctions du package cartography, plus maintenu, et qu’il est conseillé d’utiliser dorénavant le package mapsf.
Résultats de l’extraction
Dénombrement des 15 objets OSM les plus fréquents décrits par la clé amenity en Tunisie.
# Import des aménités préparées en amont
osm_pt <- st_read("data/tun/geom/tun_osm.gpkg", layer = "poi", quiet = TRUE)
# Nombre de points avec le tag "amenity".
tbl <- table(osm_pt$amenity)
tbl <- tbl[order(tbl, decreasing = TRUE)]
tbl[1:15]
cafe restaurant place_of_worship bank
2680 1272 1094 996
school pharmacy fast_food fuel
916 700 597 525
post_office police bar parking
327 153 135 135
doctors atm dentist
113 109 109 # On ne prend ici que les cliniques et hôpitaux
osm_pt <- osm_pt[osm_pt$amenity %in% c("clinic", "hospital"), ]Nombre de délégations, secteurs et cliniques/hôpitaux inclus dans OSM.
# Import des unités territoriales préparées en amont
deleg <- st_read("data/tun/geom/tun_osm.gpkg", layer = "deleg", quiet = TRUE)
sect <- st_read("data/tun/geom/tun_osm.gpkg", layer = "sect", quiet = TRUE)
# Nombre d'objets (lignes)
nrow(deleg)[1] 266nrow(sect)[1] 2104nrow(osm_pt)[1] 119Le résultat sous forme d’une carte interactive avec le package leaflet, qui sera présenté ultérieurement dans la formation.
Code
# Choisir la palette
library(leaflet)
pal <- colorFactor(palette = c("red", "gold"),
domain = c("clinic", "hospital"))
# Cartographie interactive
leaflet(osm_pt) |> # Emprise = hôpitaux
addProviderTiles("OpenStreetMap.HOT") |> # Type de tuiles chargées
addPolygons(data = sect, # Secteurs
col = "white",
fillColor = "lightgrey",
fillOpacity = 0.7,
weight = 1,
popup = paste0("<b>", sect$id_tn, "<br></b>",
sect$name_ar, "<br>", sect$name_fr),
group = "Secteurs")|>
addPolygons(data = deleg, # Délégations
col = "darkgrey",
fill = "lightgrey",
fillOpacity = 0,
weight = 1.2,
popup = paste0("<b>", deleg$id_tn, " / ",
deleg$id_hasc_deleg,"<br></b>",
deleg$name_ar, "<br>", deleg$name_fr),
group = "Délégations") |>
addCircleMarkers(radius = 4, # Hôpitaux
stroke = FALSE,
color = ~ pal(amenity),
fillOpacity = 1,
popup = paste0(osm_pt$name_ar, "<br>", osm_pt$name_fr),
group = "Cliniques et hôpitaux") |>
addLegend(pal = pal, # Légende pour différencier cliniques et hôpitaux
values = c("clinic", "hospital"),
opacity = 0.7,
title = "OSM amenity",
position = "bottomright") |>
addLayersControl(overlayGroups = c("Secteurs", "Délégations", "Cliniques et hôpitaux"),
options = layersControlOptions(collapsed = FALSE))